Survey: Light logger form factor and placement

Author

Johannes Zauner, Anna Biller, Manuel Spitschan

Preface

This is the analysis for an online survey to gather information about the acceptability of form factor and placements of light loggers. The survey was conducted using Prolific on three days in 2023, each collecting from a different sample (USA, UK, Worldwide).

This is the second part of the analysis, focusing on inference.

Import

We will start off by importing the data and setting factor levels.

Code
#Load Hmisc library
library(Hmisc)
library(tidyverse)
library(gtsummary)
library(gt)
library(ggtext)
library(rlang)
library(patchwork)
library(ordinal)
library(rlang)
library(glue)
library(broom)
library(ggridges)

#Read Data in the long format
load(file = "data/cleaned/LightLoggerFormFacto_Long_CLEANED_DATA_2024-07-02.RData")

Inferential statistics

The base hypothesis is that the ratings depend on the Wearing Position, with Sex and Sample as control variables. Further, the effect of the Wearing Position might be dependent on Sex.

The inferential analysis fits cumulative link mixed models (CLMM) to the data. The models are fitted using the ordinal package. The models are fitted for each rating separately. The base model is fitted with the following formula:

CLMM: \(E(logit(y)) = \alpha_0 + \beta_{i1} * WearingPosition + \beta_{i2} * Sex + \beta_{i3} * WearingPosition*Sex + \beta_{i4} * Sample + \alpha_j\)

Where \(E(logit(y))\) is the expected value of the rating on a logit scale, \(\alpha_0\) is the intercept, \(\beta_{i1}\) is the effect size of Wearing Position, \(\beta_{i2}\) is the effect size of Sex, \(\beta_{i3}\) is the interaction effect size of Wearing Position and Sex, \(\beta_{i4}\) is the effect size of Sample, and \(\alpha_j\) is the random intercept for each participant.

or in Wilkinson notation:

\(Rating \sim WearingPosition*Sex + Sample + (1|Id)\) (m1)

Models are compared through likelihood ratio tests. The following models are used for comparison:

A model without the interaction of Wearing Position and Sex (m2), a model without Sex (m3), a model without Sample (m4), a model with only Wearing Position (m5), and a null model, with only random effects (m0).

  • \(Rating \sim WearingPosition + Sex + Sample + (1|Id)\) (m2)

  • \(Rating \sim WearingPosition + Sample + (1|Id)\) (m3)

  • \(Rating \sim WearingPosition*Sex + (1|Id)\) (m4)

  • \(Rating \sim WearingPosition + (1|Id)\) (m5)

  • \(Rating \sim 1 + (1|Id)\) (m0)

The following comparisons are made: 1. Interaction of Sex and Wearing Position (m1 vs. m2) 2. Inclusion of Sex (m2 vs. m3) 3. Inclusion of Sample (m1 vs. m4) 4. Inclusion of Wearing Position (m5 vs. m0)

The p-values are adjusted using the Benjamini & Hochberg (1995) method, also called the false discovery rate (fdr) with an n=4 for four comparisons.

Setup

Code
#convert some variables to factor
data_long$sample_location <- as.factor(data_long$sample_location)
data_long$record_id <- as.factor(data_long$record_id)

#function generation
Infe <- function(parameter, data) {
  parameter <- enexpr(parameter)
  
  #formula generation
  formulas <- inject(c(
    #full formula
    m1 = !!parameter ~ wearing_position*demographics_sex + sample_location + (1|record_id),
    #excluding interaction from m1
    m2 = !!parameter ~ wearing_position + demographics_sex + sample_location + (1|record_id),
    #further excluding sex from m2
    m3 = !!parameter ~ wearing_position + sample_location + (1|record_id),
    #excluding sample location from m1
    m4 = !!parameter ~ wearing_position*demographics_sex + (1|record_id),
    #excluding sample location from m3
    m5 = !!parameter ~ wearing_position + (1|record_id),
    #Null model
    m0 = !!parameter ~ 1 + (1|record_id)
  ))
  
  #Model generation
  Infe_table <- tibble(
    name = names(formulas),
    formula = formulas,
    model = map(name, \(x) {switch(x,
                                        m1 = ,
                                        m2 = ,
                                        m3 = ,
                                        m4 =
                                         clmm(formulas[[x]], data = data, nAGQ = 10,
                                              subset =
                                                demographics_sex != "Other"),
                                        m5 = ,
                                        m0 = clmm(formulas[[x]], data = data, nAGQ = 10)
                                       ) 
    }
    )
  )
  
  Infe_table
}

# Function to perform anova on model pairs
anova_models <- function(pair) {
  one <- pair[[1]]
  two <- pair[[2]]
  inject(anova(Models$model[[!!one]], Models$model[[!!two]]))
}

#function to perform the model comparisons
Model_Comparisons <- function(Models) {
  # Define the model comparisons
model_comparisons <- list(
  c(1, 2),
  c(2, 3),
  c(1, 4),
  c(5, 6)
)

tibble(
  desc = c("interaction of sex and wearing position",
           "inclusion of sex",
           "inclusion of sample location",
           "inclusion of wearing position"),
  p.value = map(model_comparisons, 
                   \(x) anova_models(x) %>% .$`Pr(>Chisq)` %>% .[[2]]
                   ) %>% unlist() %>% p.adjust(method = "fdr", n = 4) %>% 
                    {case_when(
                     . < 0.001 ~ "<0.001",
                     TRUE ~ as.character(round(., 3))
                   )}
    
)
}

#probability of different answers
#setting up a function
pred <- function(eta, theta, levels, cat = 1:(length(theta)+1), inv.link = plogis) {
   Theta <- c(-Inf, theta, Inf)
   eta <- c(reference = 0, eta)
   table <- 
   sapply(cat, function(j)
          inv.link(Theta[j+1] - eta) - inv.link(Theta[j] - eta) ) 
   colnames(table) <- levels
   table
}

data_significance_matrix <- function(Model, subset = "none") {
  Models <- tibble(
  reference = levels(data_long$wearing_position),
  data = map(reference, 
             \(x) data_long %>% mutate(wearing_position = 
                                         fct_relevel(wearing_position, x))
             ),
  model = map(data, 
              \(x) switch(subset,
                          none = clmm(Model$formula, data = x, nAGQ = 10),
                          clmm(Model$formula, data = x, nAGQ = 10,
                               subset = demographics_sex != "Other")
                          )
              ),
  tidy = map(model, 
             \(x) tidy(x) %>% 
               filter(str_detect(term, "wearing_position")) %>% 
               select(term, p.value) %>% 
               mutate(term = str_remove(term, "wearing_position"))
             )
  )

Models %>% 
  select(-data, -model) %>% 
  unnest(tidy) %>% 
  rowwise() %>% 
  mutate(
    p.value = p.adjust(p.value, method = "fdr", n = 7),
    different = p.value <= 0.05) %>%
  ungroup() %>% 
  rbind(
    tibble(reference = levels(data_long$wearing_position),
           term = reference,
           p.value = NA,
           different = FALSE))
}

#McFadden´s Pseudo R2, values between 0.2 bis 0.4 are an excellent model fit
McFadden <- function(model, nullmodel) {
  ll.null <- nullmodel$logLik
  ll.proposed <- model$logLik
  1 - ll.proposed/ll.null
}

#Formula to plot the significance matrix
Significance_matrix <- function(sig_matrix, sex = FALSE){
  sig_matrix %>% 
  mutate(
    across(c(reference, term), \(x) fct_relevel(x, levels(data_long$wearing_position)))) %>%
  ggplot(aes(x=reference, y = term, fill = different)) +
  geom_tile(col = "grey30") + 
  scale_fill_manual(values = c("white", "skyblue4"))+
  coord_fixed() +
  theme_void() +
  labs(title = label(data_long[parameter]),
       subtitle = case_when(!sex ~ "<span style='color:skyblue4'>Significant differences</span> due to wearing position",
                            sex ~ "Significant differences due to <span style='color:skyblue4'>wearing position</span> and <span style='color:red'>sex</span>" )
       )+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text.y = element_text(hjust = 1),
        plot.title = 
          element_textbox_simple(
            hjust = 0,
            padding = margin(0.5, 0, 0.5, 0, "cm")),
        plot.subtitle = element_textbox_simple(),
        plot.title.position = "plot",
        )
}

#table of probabilities
prob_table <- function(beta, title){
pred(beta, Model$Theta, 
     levels = levels(data_long[[parameter]])) %>% 
  as_tibble(rownames = "Wearing Position") %>%
  mutate(`Wearing Position` = 
           str_remove(`Wearing Position`, "wearing_position"),
         `Wearing Position` = str_replace(`Wearing Position`, "reference", "Chest Pin")
           ) %>%
  gt() %>% fmt_percent() %>% 
  gt::tab_caption(md(glue(title)))
}

Model selection

Code
parameter <- "appearance"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for appearance
desc p.value
interaction of sex and wearing position 0.827
inclusion of sex 0.827
inclusion of sample location 0.57
inclusion of wearing position <0.001

The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: appearance ~ wearing_position + (1 | record_id)
data:    data

 link  threshold nobs logLik   AIC     niter      max.grad cond.H 
 logit flexible  1160 -1774.52 3577.05 1604(7930) 2.38e-03 1.6e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 1.531    1.237   
Number of groups:  record_id 145 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
wearing_positionWrist          -0.4594     0.2146  -2.141   0.0323 *  
wearing_positionNecklace       -1.0217     0.2161  -4.728 2.26e-06 ***
wearing_positionSleeve collar  -1.5935     0.2253  -7.072 1.53e-12 ***
wearing_positionCollar pin     -2.2585     0.2212 -10.211  < 2e-16 ***
wearing_positionNeck loop      -2.2304     0.2249  -9.919  < 2e-16 ***
wearing_positionHat pin        -2.5887     0.2239 -11.563  < 2e-16 ***
wearing_positionGlasses        -3.1912     0.2349 -13.584  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                         Estimate Std. Error z value
Extremely unattractive|Very unattractive  -4.6499     0.2330 -19.957
Very unattractive|Unattractive            -3.7687     0.2179 -17.295
Unattractive|Neutral                      -2.3197     0.2008 -11.553
Neutral|Attractive                        -0.6844     0.1886  -3.629
Attractive|Very attractive                 1.7319     0.2029   8.537
Very attractive|Extremely attractive       3.3479     0.2744  12.202
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.08141085
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model)

P1 <- Significance_matrix(sig_matrix)

ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P1, width = 4, height = 4.6)

P1

Predicted probabilities

Code
#if one wanted to print predictions for different subject percentiles
stDev <-  as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev, 
           "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for appearance
Wearing Position Extremely unattractive Very unattractive Unattractive Neutral Attractive Very attractive Extremely attractive
Chest Pin 0.95% 1.31% 6.69% 24.58% 51.44% 11.64% 3.40%
Wrist 0.30% 0.42% 2.28% 10.69% 50.30% 25.95% 10.05%
Necklace 0.07% 0.10% 0.57% 2.96% 26.44% 38.33% 31.52%
Sleeve collar 0.02% 0.02% 0.14% 0.73% 8.38% 24.71% 66.01%
Collar pin 0.00% 0.00% 0.03% 0.14% 1.71% 6.92% 91.20%
Neck loop 0.00% 0.00% 0.03% 0.15% 1.83% 7.37% 90.61%
Hat pin 0.00% 0.00% 0.01% 0.06% 0.75% 3.21% 95.97%
Glasses 0.00% 0.00% 0.00% 0.01% 0.17% 0.73% 99.09%
Code
#average:
prob_table(Model$beta, 
           "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for appearance
Wearing Position Extremely unattractive Very unattractive Unattractive Neutral Attractive Very attractive Extremely attractive
Chest Pin 0.95% 1.31% 6.69% 24.58% 51.44% 11.64% 3.40%
Wrist 1.49% 2.03% 9.94% 30.93% 45.55% 7.88% 2.17%
Necklace 2.59% 3.44% 15.42% 36.90% 35.66% 4.74% 1.25%
Sleeve collar 4.49% 5.71% 22.40% 38.68% 25.25% 2.76% 0.71%
Collar pin 8.38% 9.71% 30.38% 34.37% 15.35% 1.45% 0.37%
Neck loop 8.17% 9.51% 30.09% 34.66% 15.70% 1.49% 0.38%
Hat pin 11.29% 12.21% 33.18% 30.35% 11.65% 1.05% 0.26%
Glasses 18.87% 17.08% 34.55% 21.96% 6.82% 0.58% 0.14%
Code
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev, 
           "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for appearance
Wearing Position Extremely unattractive Very unattractive Unattractive Neutral Attractive Very attractive Extremely attractive
Chest Pin 0.95% 1.31% 6.69% 24.58% 51.44% 11.64% 3.40%
Wrist 2.95% 3.89% 16.97% 37.78% 33.14% 4.18% 1.09%
Necklace 11.13% 12.08% 33.07% 30.57% 11.82% 1.07% 0.27%
Sleeve collar 34.56% 21.48% 28.40% 12.09% 3.15% 0.26% 0.06%
Collar pin 73.80% 13.38% 9.48% 2.67% 0.61% 0.05% 0.01%
Neck loop 72.41% 13.96% 10.06% 2.86% 0.65% 0.05% 0.01%
Hat pin 86.61% 7.37% 4.54% 1.19% 0.27% 0.02% 0.01%
Glasses 96.72% 1.89% 1.06% 0.26% 0.06% 0.00% 0.00%

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)

This sections analyses how the context of wear influences the ratings, dependent on the wearing position

Code
#function to remove labels. This is necessary for pivot_longer to work
clear.labels <- function(x) {
  if(is.list(x)) {
    for(i in seq_along(x)) {
      class(x[[i]]) <- setdiff(class(x[[i]]), 'labelled') 
      attr(x[[i]],"label") <- NULL
    } 
  } else {
    class(x) <- setdiff(class(x), "labelled")
    attr(x, "label") <- NULL
  }
  return(x)
}

#pivoting the data even longer to get the context into the observation
data_long2 <- data_long %>% 
  select(record_id,
         demographics_sex,
         wearing_position,
         starts_with("wearlikelihoodrating")
  ) %>% 
  clear.labels() %>% 
  pivot_longer(-c(record_id, demographics_sex, wearing_position), 
               names_to = "context", 
               values_to = "rating") %>% 
  mutate(context = 
           str_remove(context, "wearlikelihoodrating_") %>% 
           factor() %>% 
           fct_reorder(as.numeric(rating), .desc = TRUE, .fun = mean)
           # fct_relevel(c("public", "work", "home", "social", "exercise")
           )

#model generation
Model <- clmm(rating ~ context*wearing_position + (1|record_id), data = data_long2, nAGQ = 10)
Model2 <- clmm(rating ~ context+wearing_position + (1|record_id), data = data_long2, nAGQ = 10)
Model3 <- clmm(rating ~ context + (1|record_id), data = data_long2, nAGQ = 10)
Model4 <- clmm(rating ~ wearing_position + (1|record_id), data = data_long2, nAGQ = 10)
Model0 <- clmm(rating ~ 1 + (1|record_id), data = data_long2, nAGQ = 10)

#model selection
anova(Model, Model2)
Likelihood ratio tests of cumulative link models:
 
       formula:                                              link: threshold:
Model2 rating ~ context + wearing_position + (1 | record_id) logit flexible  
Model  rating ~ context * wearing_position + (1 | record_id) logit flexible  

       no.par   AIC  logLik LR.stat df Pr(>Chisq)    
Model2     18 18590 -9276.8                          
Model      46 18448 -9178.1  197.41 28  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
anova(Model2, Model3)
Likelihood ratio tests of cumulative link models:
 
       formula:                                              link: threshold:
Model3 rating ~ context + (1 | record_id)                    logit flexible  
Model2 rating ~ context + wearing_position + (1 | record_id) logit flexible  

       no.par   AIC  logLik LR.stat df Pr(>Chisq)    
Model3     11 19937 -9957.5                          
Model2     18 18590 -9276.8  1361.5  7  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
anova(Model3, Model0)
Likelihood ratio tests of cumulative link models:
 
       formula:                           link: threshold:
Model0 rating ~ 1 + (1 | record_id)       logit flexible  
Model3 rating ~ context + (1 | record_id) logit flexible  

       no.par   AIC   logLik LR.stat df Pr(>Chisq)    
Model0      7 20120 -10053.2                          
Model3     11 19937  -9957.5  191.28  4  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
anova(Model4, Model0)
Likelihood ratio tests of cumulative link models:
 
       formula:                                    link: threshold:
Model0 rating ~ 1 + (1 | record_id)                logit flexible  
Model4 rating ~ wearing_position + (1 | record_id) logit flexible  

       no.par   AIC   logLik LR.stat df Pr(>Chisq)    
Model0      7 20120 -10053.2                          
Model4     14 18831  -9401.6  1303.1  7  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
#model summaries
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: rating ~ context * wearing_position + (1 | record_id)
data:    data_long2

 link  threshold nobs logLik   AIC      niter        max.grad cond.H 
 logit flexible  5800 -9178.09 18448.18 12392(93749) 6.64e-03 3.9e+03

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 2.026    1.423   
Number of groups:  record_id 145 

Coefficients:
                                              Estimate Std. Error z value
contextpublic                                  -0.4639     0.2067  -2.244
contextexercise                                -0.5336     0.2125  -2.511
contextwork                                    -0.7212     0.2086  -3.458
contextsocial                                  -0.9508     0.2078  -4.576
wearing_positionWrist                          -0.7158     0.2148  -3.333
wearing_positionNecklace                       -0.6401     0.2137  -2.996
wearing_positionSleeve collar                  -1.4403     0.2170  -6.636
wearing_positionCollar pin                     -1.3705     0.2114  -6.481
wearing_positionNeck loop                      -1.7915     0.2148  -8.338
wearing_positionHat pin                        -2.9718     0.2167 -13.713
wearing_positionGlasses                        -2.4774     0.2142 -11.569
contextpublic:wearing_positionWrist             0.3498     0.2977   1.175
contextexercise:wearing_positionWrist           0.8694     0.3055   2.846
contextwork:wearing_positionWrist               0.2889     0.3021   0.956
contextsocial:wearing_positionWrist             0.3082     0.2979   1.035
contextpublic:wearing_positionNecklace         -0.1876     0.2965  -0.633
contextexercise:wearing_positionNecklace       -1.3970     0.3048  -4.584
contextwork:wearing_positionNecklace           -0.3569     0.2996  -1.191
contextsocial:wearing_positionNecklace         -0.2508     0.2978  -0.842
contextpublic:wearing_positionSleeve collar    -0.2226     0.3019  -0.737
contextexercise:wearing_positionSleeve collar   1.3441     0.3107   4.326
contextwork:wearing_positionSleeve collar      -0.5646     0.3007  -1.877
contextsocial:wearing_positionSleeve collar    -0.4524     0.3004  -1.506
contextpublic:wearing_positionCollar pin       -0.5273     0.2931  -1.799
contextexercise:wearing_positionCollar pin      0.1147     0.3011   0.381
contextwork:wearing_positionCollar pin         -0.6189     0.2981  -2.076
contextsocial:wearing_positionCollar pin       -0.6357     0.2964  -2.145
contextpublic:wearing_positionNeck loop        -0.2973     0.2964  -1.003
contextexercise:wearing_positionNeck loop      -0.3186     0.3064  -1.040
contextwork:wearing_positionNeck loop          -0.5299     0.3005  -1.763
contextsocial:wearing_positionNeck loop        -0.3003     0.2990  -1.005
contextpublic:wearing_positionHat pin           0.5526     0.2968   1.862
contextexercise:wearing_positionHat pin         0.3640     0.3055   1.192
contextwork:wearing_positionHat pin            -0.1569     0.3023  -0.519
contextsocial:wearing_positionHat pin           0.4544     0.3003   1.513
contextpublic:wearing_positionGlasses          -0.2479     0.2995  -0.828
contextexercise:wearing_positionGlasses        -0.7585     0.3049  -2.488
contextwork:wearing_positionGlasses            -0.1357     0.3020  -0.449
contextsocial:wearing_positionGlasses          -0.2515     0.3036  -0.829
                                              Pr(>|z|)    
contextpublic                                 0.024811 *  
contextexercise                               0.012028 *  
contextwork                                   0.000545 ***
contextsocial                                 4.74e-06 ***
wearing_positionWrist                         0.000860 ***
wearing_positionNecklace                      0.002735 ** 
wearing_positionSleeve collar                 3.21e-11 ***
wearing_positionCollar pin                    9.09e-11 ***
wearing_positionNeck loop                      < 2e-16 ***
wearing_positionHat pin                        < 2e-16 ***
wearing_positionGlasses                        < 2e-16 ***
contextpublic:wearing_positionWrist           0.239885    
contextexercise:wearing_positionWrist         0.004433 ** 
contextwork:wearing_positionWrist             0.338828    
contextsocial:wearing_positionWrist           0.300884    
contextpublic:wearing_positionNecklace        0.526902    
contextexercise:wearing_positionNecklace      4.56e-06 ***
contextwork:wearing_positionNecklace          0.233557    
contextsocial:wearing_positionNecklace        0.399794    
contextpublic:wearing_positionSleeve collar   0.460819    
contextexercise:wearing_positionSleeve collar 1.52e-05 ***
contextwork:wearing_positionSleeve collar     0.060497 .  
contextsocial:wearing_positionSleeve collar   0.132102    
contextpublic:wearing_positionCollar pin      0.072055 .  
contextexercise:wearing_positionCollar pin    0.703195    
contextwork:wearing_positionCollar pin        0.037868 *  
contextsocial:wearing_positionCollar pin      0.031984 *  
contextpublic:wearing_positionNeck loop       0.315737    
contextexercise:wearing_positionNeck loop     0.298500    
contextwork:wearing_positionNeck loop         0.077833 .  
contextsocial:wearing_positionNeck loop       0.315105    
contextpublic:wearing_positionHat pin         0.062620 .  
contextexercise:wearing_positionHat pin       0.233364    
contextwork:wearing_positionHat pin           0.603725    
contextsocial:wearing_positionHat pin         0.130261    
contextpublic:wearing_positionGlasses         0.407884    
contextexercise:wearing_positionGlasses       0.012838 *  
contextwork:wearing_positionGlasses           0.653081    
contextsocial:wearing_positionGlasses         0.407351    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -4.2432     0.1984 -21.386
Very Unlikely|Unlikely            -3.3268     0.1960 -16.973
Unlikely|Neutral                  -2.2250     0.1940 -11.470
Neutral|Likely                    -1.2738     0.1927  -6.609
Likely|Very likely                 0.4591     0.1917   2.395
Very likely|Extremely likely       2.1776     0.1976  11.021
Code
summary(Model3)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: rating ~ context + (1 | record_id)
data:    data_long2

 link  threshold nobs logLik   AIC      niter      max.grad cond.H 
 logit flexible  5800 -9957.53 19937.05 1209(8342) 8.62e-03 3.6e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 1.451    1.205   
Number of groups:  record_id 145 

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
contextpublic   -0.45146    0.07412  -6.091 1.12e-09 ***
contextexercise -0.44588    0.07606  -5.862 4.57e-09 ***
contextwork     -0.83751    0.07543 -11.103  < 2e-16 ***
contextsocial   -0.91921    0.07517 -12.229  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -2.3870     0.1188 -20.100
Very Unlikely|Unlikely            -1.5906     0.1163 -13.674
Unlikely|Neutral                  -0.6642     0.1149  -5.782
Neutral|Likely                     0.1109     0.1145   0.968
Likely|Very likely                 1.5633     0.1167  13.391
Very likely|Extremely likely       3.1163     0.1288  24.187
Code
summary(Model4)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: rating ~ wearing_position + (1 | record_id)
data:    data_long2

 link  threshold nobs logLik   AIC      niter       max.grad cond.H 
 logit flexible  5800 -9401.64 18831.28 2076(16279) 9.58e-03 4.0e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 1.851    1.36    
Number of groups:  record_id 145 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
wearing_positionWrist         -0.33834    0.09403  -3.598 0.000321 ***
wearing_positionNecklace      -1.01892    0.09474 -10.755  < 2e-16 ***
wearing_positionSleeve collar -1.39982    0.09718 -14.405  < 2e-16 ***
wearing_positionCollar pin    -1.62095    0.09572 -16.935  < 2e-16 ***
wearing_positionNeck loop     -1.99074    0.09701 -20.522  < 2e-16 ***
wearing_positionHat pin       -2.59863    0.09921 -26.194  < 2e-16 ***
wearing_positionGlasses       -2.62203    0.10032 -26.136  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -3.5293     0.1392 -25.352
Very Unlikely|Unlikely            -2.6479     0.1362 -19.441
Unlikely|Neutral                  -1.6022     0.1338 -11.972
Neutral|Likely                    -0.7050     0.1324  -5.323
Likely|Very likely                 0.9436     0.1327   7.109
Very likely|Extremely likely       2.6047     0.1431  18.196
Code
#get the predicted probabilities depending on context, save as a table
context_preds <- pred(Model3$beta, Model3$Theta, 
     levels = levels(data_long2$rating)) %>% 
  as_tibble(rownames = "Parameters") %>%
  mutate(`Parameters` = 
           str_remove(`Parameters`, "wearing_position|context"),
         `Parameters` = str_replace(`Parameters`, "reference", "home")
           )

context_preds_table <- 
context_preds %>% 
  gt() %>% 
  fmt_percent(decimals = 0) %>% 
  gt::tab_caption(md("Likelihood of wear ratings depending on context")) %>% 
  gt::data_color(columns = -Parameters, method = "numeric", 
                 colors = scales::col_numeric(palette = "viridis", 
                                              domain = c(0,0.30))
                 ) %>% 
  cols_width(everything() ~ px(100))
Warning: Since gt v0.9.0, the `colors` argument has been deprecated.
• Please use the `fn` argument instead.
This warning is displayed once every 8 hours.
Code
context_preds_table
Likelihood of wear ratings depending on context
Parameters Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
home 8% 9% 17% 19% 30% 13% 4%
public 13% 12% 20% 19% 25% 9% 3%
exercise 13% 12% 20% 19% 25% 9% 3%
work 18% 14% 22% 18% 20% 6% 2%
social 19% 15% 23% 17% 19% 6% 2%
Code
gtsave(context_preds_table, "output/02_tables/table_2_context.png")
file:////var/folders/9p/326_k3kx43qbn_cyl1rqfhb00000gn/T//Rtmp2v7SBD/file85a86024d65a.html screenshot completed
Code
#get the predicted probabilities depending on wearing_position, save as a table
context_preds <- pred(Model4$beta, Model4$Theta, 
     levels = levels(data_long2$rating)) %>% 
  as_tibble(rownames = "Parameters") %>%
  mutate(`Parameters` = 
           str_remove(`Parameters`, "wearing_position|context"),
         `Parameters` = str_replace(`Parameters`, "reference", "Chest pin")
           )

context_preds_table <- 
context_preds %>% 
  gt() %>% 
  fmt_percent(decimals = 0) %>% 
  gt::tab_caption(md("Likelihood of wear ratings depending on wearing position")) %>% 
  gt::data_color(columns = -Parameters, method = "numeric", 
                 colors = scales::col_numeric(palette = "viridis", 
                                              domain = c(0,0.39))
                 ) %>% 
  cols_width(everything() ~ px(100))

context_preds_table
Likelihood of wear ratings depending on wearing position
Parameters Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest pin 3% 4% 10% 16% 39% 21% 7%
Wrist 4% 5% 13% 19% 37% 17% 5%
Necklace 8% 9% 19% 22% 30% 10% 3%
Sleeve collar 11% 12% 23% 22% 25% 7% 2%
Collar pin 13% 13% 24% 21% 21% 6% 1%
Neck loop 18% 16% 25% 19% 17% 4% 1%
Hat pin 28% 20% 24% 14% 10% 2% 1%
Glasses 29% 21% 24% 14% 10% 2% 1%
Code
gtsave(context_preds_table, "output/02_tables/table_3_wearing_position.png")
file:////var/folders/9p/326_k3kx43qbn_cyl1rqfhb00000gn/T//Rtmp2v7SBD/file85a86e794f42.html screenshot completed
Code
#create a table of the interaction model
#first, collect all the beta coefficients for the context
coefs_context <- 
  Model$beta %>% 
  as_tibble(rownames = "Parameters") %>% 
  filter(str_detect(Parameters, "context"),
         !str_detect(Parameters, "wearing_position")) %>% 
  rename(Context = Parameters,
         beta1 = value) %>% 
  add_row(Context = "contexthome", beta1 = 0, .before = 1)

#second, collect all the beta coefficients for the wearing position
coefs_position <- 
  Model$beta %>% 
  as_tibble(rownames = "Parameters") %>% 
  filter(str_detect(Parameters, "wearing_position"),
         !str_detect(Parameters, "context")) %>% 
    rename(Position2 = Parameters,
         beta2 = value) %>% 
  add_row(Position2 = "wearing_positionChest pin", beta2 = 0, .before = 1)

#third, collect all the beta coefficients for the interaction
coefs_extra <- 
  Model$beta %>% 
  as_tibble(rownames = "Parameters") %>% 
  filter(str_detect(Parameters, "wearing_position")) %>% 
  rename(Position = Parameters,
         beta3 = value) %>% 
  add_row(Position = "contextpublic:wearing_positionChest pin", beta3 = 0, .before = 2) %>%
  add_row(Position = "contextexercise:wearing_positionChest pin", beta3 = 0, .before = 3) %>%
  add_row(Position = "contextwork:wearing_positionChest pin", beta3 = 0, .before = 4) %>%
  add_row(Position = "contextsocial:wearing_positionChest pin", beta3 = 0, .before = 5) %>%
  mutate(Position = 
           str_replace(
             Position, "^wearing_position", "contexthome:wearing_position")
         )

#combine the three tables to a single table and filter out non-relevant combinations
coefs_combined <- 
expand_grid(coefs_context, coefs_position, coefs_extra) %>% 
  filter(str_detect(Position, Context),
         str_detect(Position, Position2)) %>% 
  mutate(Context = str_remove(Context, "context"), 
         Position = str_remove(Position, "(.*):wearing_position"),
         beta3 = ifelse(str_detect(Context , "home"), 0, beta3))

coefs_combined <- 
  set_names(coefs_combined$beta1+coefs_combined$beta2 + coefs_combined$beta3, 
            paste(coefs_combined$Context, coefs_combined$Position, sep = ":"))

#create and save a table of the predicted rating
combined_pred <- 
pred(coefs_combined, Model$Theta, 
     levels = levels(data_long2$rating)
     ) %>% 
  as_tibble(rownames = "Parameters") %>%
  mutate(
         `Parameters` = str_replace(`Parameters`, "reference", "home:Chest pin")
  ) %>% 
  separate_wider_delim(Parameters, ":", names = c("Context", "Position")) %>% 
  pivot_longer(cols = -c(Context, Position), names_to = "Rating", values_to = "Probability") %>% 
  mutate(Rating = factor(Rating, levels = levels(data_long2$rating))) %>% 
  group_by(Context, Position) %>% 
  summarize(Probability = Rating[which.max(Probability)], .groups = "drop") %>% 
  pivot_wider(names_from = Position, values_from = Probability) %>% 
  arrange(factor(Context, levels = levels(data_long2$context))) %>% 
  select(Context, all_of(levels(data_long2$wearing_position))) %>% 
  gt() %>% 
  gt::tab_caption(md("*Likelihood of wear* rating prediction depending on **context** and **wearing position**")) %>% 
  gt::data_color(columns = -Context, ordered = TRUE,
                 palette = "viridis"
                 ) %>% 
  cols_width(everything() ~ px(100))

combined_pred
Likelihood of wear rating prediction depending on context and wearing position
Context Chest pin Wrist Necklace Sleeve collar Collar pin Neck loop Hat pin Glasses
home Likely Likely Likely Likely Likely Likely Unlikely Unlikely
public Likely Likely Likely Unlikely Unlikely Unlikely Unlikely Extremely unlikely
exercise Likely Likely Unlikely Likely Likely Unlikely Unlikely Extremely unlikely
work Likely Likely Likely Unlikely Unlikely Unlikely Extremely unlikely Extremely unlikely
social Likely Likely Likely Unlikely Unlikely Unlikely Extremely unlikely Extremely unlikely
Code
gtsave(combined_pred, "output/02_tables/table_4_context_position.png")
file:////var/folders/9p/326_k3kx43qbn_cyl1rqfhb00000gn/T//Rtmp2v7SBD/file85a85e42aeb7.html screenshot completed

Model selection

Code
parameter <- "wearlikelihoodrating_public"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for wearlikelihoodrating_public
desc p.value
interaction of sex and wearing position 0.812
inclusion of sex 0.812
inclusion of sample location 0.389
inclusion of wearing position <0.001

The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: wearlikelihoodrating_public ~ wearing_position + (1 | record_id)
data:    data

 link  threshold nobs logLik   AIC     niter      max.grad cond.H 
 logit flexible  1160 -1899.19 3826.39 1570(9364) 1.01e-03 2.2e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 1.921    1.386   
Number of groups:  record_id 145 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
wearing_positionWrist          -0.3931     0.2132  -1.844   0.0652 .  
wearing_positionNecklace       -0.8693     0.2138  -4.067 4.77e-05 ***
wearing_positionSleeve collar  -1.7525     0.2218  -7.901 2.76e-15 ***
wearing_positionCollar pin     -2.0125     0.2176  -9.250  < 2e-16 ***
wearing_positionNeck loop      -2.1992     0.2195 -10.021  < 2e-16 ***
wearing_positionHat pin        -2.5519     0.2234 -11.425  < 2e-16 ***
wearing_positionGlasses        -2.8638     0.2317 -12.358  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -4.0037     0.2262 -17.699
Very Unlikely|Unlikely            -3.1317     0.2140 -14.633
Unlikely|Neutral                  -1.9030     0.2022  -9.410
Neutral|Likely                    -0.9430     0.1958  -4.816
Likely|Very likely                 1.0424     0.1964   5.307
Very likely|Extremely likely       2.8424     0.2398  11.855
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.07070644
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model)

P2 <- Significance_matrix(sig_matrix)

ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P2, width = 4, height = 4.7)

P2

Predicted probabilities

Code
#if one wanted to print predictions for different subject percentiles
stDev <-  as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev, 
           "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_public
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 1.79% 2.39% 8.79% 15.05% 45.90% 20.56% 5.51%
Wrist 0.52% 0.72% 2.88% 5.98% 34.91% 38.18% 16.80%
Necklace 0.12% 0.16% 0.67% 1.49% 12.95% 36.99% 47.63%
Sleeve collar 0.01% 0.01% 0.04% 0.09% 0.95% 5.22% 93.68%
Collar pin 0.00% 0.00% 0.02% 0.04% 0.42% 2.39% 97.12%
Neck loop 0.00% 0.00% 0.01% 0.02% 0.23% 1.35% 98.38%
Hat pin 0.00% 0.00% 0.00% 0.01% 0.08% 0.45% 99.46%
Glasses 0.00% 0.00% 0.00% 0.00% 0.03% 0.17% 99.80%
Code
#average:
prob_table(Model$beta, 
           "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_public
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 1.79% 2.39% 8.79% 15.05% 45.90% 20.56% 5.51%
Wrist 2.63% 3.44% 12.02% 18.49% 44.19% 15.44% 3.79%
Necklace 4.17% 5.26% 16.81% 21.92% 38.96% 10.49% 2.39%
Sleeve collar 9.52% 10.59% 26.13% 22.96% 25.04% 4.76% 1.00%
Collar pin 12.01% 12.60% 28.12% 21.72% 21.05% 3.73% 0.77%
Neck loop 14.13% 14.11% 29.11% 20.49% 18.40% 3.12% 0.64%
Hat pin 18.97% 16.93% 29.78% 17.65% 14.00% 2.22% 0.45%
Glasses 24.23% 19.11% 28.99% 14.89% 10.80% 1.64% 0.33%
Code
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev, 
           "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_public
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 1.79% 2.39% 8.79% 15.05% 45.90% 20.56% 5.51%
Wrist 5.94% 7.19% 20.93% 23.37% 33.33% 7.59% 1.66%
Necklace 22.16% 18.35% 29.43% 15.93% 11.92% 1.84% 0.37%
Sleeve collar 82.27% 9.46% 5.70% 1.57% 0.86% 0.12% 0.02%
Collar pin 91.35% 4.84% 2.66% 0.70% 0.38% 0.05% 0.01%
Neck loop 95.01% 2.84% 1.51% 0.39% 0.21% 0.03% 0.01%
Hat pin 98.31% 0.98% 0.50% 0.13% 0.07% 0.01% 0.00%
Glasses 99.36% 0.37% 0.19% 0.05% 0.03% 0.00% 0.00%

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)

Model selection

Code
parameter <- "wearlikelihoodrating_work"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for wearlikelihoodrating_work
desc p.value
interaction of sex and wearing position 0.118
inclusion of sex 0.648
inclusion of sample location 0.077
inclusion of wearing position <0.001

The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: wearlikelihoodrating_work ~ wearing_position + (1 | record_id)
data:    data

 link  threshold nobs logLik   AIC     niter      max.grad cond.H 
 logit flexible  1160 -1874.88 3777.76 1567(7749) 1.18e-03 2.1e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 2.229    1.493   
Number of groups:  record_id 145 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
wearing_positionWrist          -0.4633     0.2148  -2.157    0.031 *  
wearing_positionNecklace       -1.0287     0.2140  -4.806 1.54e-06 ***
wearing_positionSleeve collar  -2.0777     0.2201  -9.441  < 2e-16 ***
wearing_positionCollar pin     -2.0459     0.2188  -9.351  < 2e-16 ***
wearing_positionNeck loop      -2.4086     0.2227 -10.818  < 2e-16 ***
wearing_positionHat pin        -3.2622     0.2331 -13.997  < 2e-16 ***
wearing_positionGlasses        -2.7678     0.2303 -12.019  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -3.7448     0.2260 -16.570
Very Unlikely|Unlikely            -2.6696     0.2135 -12.506
Unlikely|Neutral                  -1.5506     0.2046  -7.577
Neutral|Likely                    -0.5348     0.1989  -2.688
Likely|Very likely                 1.2863     0.2036   6.317
Very likely|Extremely likely       2.9820     0.2503  11.913
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.08277394
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model)

P3 <- Significance_matrix(sig_matrix)

ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P3, width = 4, height = 4.7)

P3

Predicted probabilities

Code
#if one wanted to print predictions for different subject percentiles
stDev <-  as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev, 
           "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_work
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 2.31% 4.17% 11.02% 19.44% 41.41% 16.82% 4.82%
Wrist 0.43% 0.82% 2.48% 5.94% 30.15% 38.47% 21.70%
Necklace 0.05% 0.10% 0.33% 0.84% 6.36% 23.53% 68.78%
Sleeve collar 0.00% 0.00% 0.01% 0.02% 0.15% 0.78% 99.04%
Collar pin 0.00% 0.00% 0.01% 0.02% 0.17% 0.88% 98.92%
Neck loop 0.00% 0.00% 0.00% 0.01% 0.04% 0.23% 99.71%
Hat pin 0.00% 0.00% 0.00% 0.00% 0.00% 0.01% 99.99%
Glasses 0.00% 0.00% 0.00% 0.00% 0.01% 0.06% 99.92%
Code
#average:
prob_table(Model$beta, 
           "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_work
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 2.31% 4.17% 11.02% 19.44% 41.41% 16.82% 4.82%
Wrist 3.62% 6.30% 15.29% 23.00% 36.98% 11.72% 3.09%
Necklace 6.20% 10.03% 21.01% 24.86% 28.91% 7.21% 1.78%
Sleeve collar 15.88% 19.74% 27.26% 19.51% 14.27% 2.71% 0.63%
Collar pin 15.46% 19.43% 27.24% 19.79% 14.63% 2.80% 0.65%
Neck loop 20.81% 22.70% 26.71% 16.46% 10.89% 1.97% 0.45%
Hat pin 38.16% 26.23% 20.31% 9.16% 5.09% 0.85% 0.19%
Glasses 27.35% 25.11% 24.70% 13.16% 7.98% 1.39% 0.32%
Code
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev, 
           "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_work
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 2.31% 4.17% 11.02% 19.44% 41.41% 16.82% 4.82%
Wrist 11.45% 16.03% 26.23% 22.51% 18.98% 3.89% 0.92%
Necklace 50.68% 24.39% 15.14% 6.00% 3.15% 0.52% 0.12%
Sleeve collar 97.96% 1.33% 0.47% 0.15% 0.07% 0.01% 0.00%
Collar pin 97.72% 1.49% 0.53% 0.17% 0.08% 0.01% 0.00%
Neck loop 99.39% 0.40% 0.14% 0.04% 0.02% 0.00% 0.00%
Hat pin 99.97% 0.02% 0.01% 0.00% 0.00% 0.00% 0.00%
Glasses 99.83% 0.11% 0.04% 0.01% 0.01% 0.00% 0.00%

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)

Model selection

Code
parameter <- "wearlikelihoodrating_home"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for wearlikelihoodrating_home
desc p.value
interaction of sex and wearing position 0.332
inclusion of sex 0.653
inclusion of sample location 0.943
inclusion of wearing position <0.001

The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: wearlikelihoodrating_home ~ wearing_position + (1 | record_id)
data:    data

 link  threshold nobs logLik   AIC     niter      max.grad cond.H 
 logit flexible  1160 -1861.86 3751.72 1483(8824) 2.94e-03 2.0e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 3.241    1.8     
Number of groups:  record_id 145 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
wearing_positionWrist          -0.7833     0.2146  -3.650 0.000263 ***
wearing_positionNecklace       -0.6834     0.2138  -3.197 0.001390 ** 
wearing_positionSleeve collar  -1.5234     0.2205  -6.908 4.92e-12 ***
wearing_positionCollar pin     -1.4878     0.2170  -6.857 7.04e-12 ***
wearing_positionNeck loop      -1.9235     0.2213  -8.691  < 2e-16 ***
wearing_positionHat pin        -3.1911     0.2339 -13.645  < 2e-16 ***
wearing_positionGlasses        -2.6719     0.2275 -11.746  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -4.4846     0.2543 -17.634
Very Unlikely|Unlikely            -3.7356     0.2432 -15.358
Unlikely|Neutral                  -2.5400     0.2303 -11.028
Neutral|Likely                    -1.4360     0.2225  -6.454
Likely|Very likely                 0.4905     0.2180   2.250
Very likely|Extremely likely       2.3480     0.2383   9.855
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.07320812
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model)

P4 <- Significance_matrix(sig_matrix)

ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P4, width = 4, height = 4.7)

P4

Predicted probabilities

Code
#if one wanted to print predictions for different subject percentiles
stDev <-  as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev, 
           "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_home
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 1.12% 1.21% 4.98% 11.91% 42.81% 29.26% 8.72%
Wrist 0.02% 0.02% 0.08% 0.24% 2.08% 11.40% 86.16%
Necklace 0.03% 0.03% 0.14% 0.41% 3.48% 17.39% 78.51%
Sleeve collar 0.00% 0.00% 0.00% 0.00% 0.04% 0.26% 99.69%
Collar pin 0.00% 0.00% 0.00% 0.01% 0.05% 0.32% 99.63%
Neck loop 0.00% 0.00% 0.00% 0.00% 0.00% 0.03% 99.96%
Hat pin 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 100.00%
Glasses 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 100.00%
Code
#average:
prob_table(Model$beta, 
           "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_home
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 1.12% 1.21% 4.98% 11.91% 42.81% 29.26% 8.72%
Wrist 2.41% 2.55% 9.76% 19.52% 43.90% 17.68% 4.18%
Necklace 2.19% 2.33% 9.00% 18.52% 44.36% 19.01% 4.60%
Sleeve collar 4.92% 4.95% 16.70% 25.62% 36.04% 9.73% 2.04%
Collar pin 4.76% 4.80% 16.33% 25.42% 36.55% 10.04% 2.11%
Neck loop 7.17% 6.87% 21.02% 26.89% 29.84% 6.83% 1.38%
Hat pin 21.53% 15.19% 29.01% 19.53% 12.28% 2.06% 0.39%
Glasses 14.03% 11.63% 27.63% 24.19% 18.45% 3.40% 0.66%
Code
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev, 
           "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_home
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 1.12% 1.21% 4.98% 11.91% 42.81% 29.26% 8.72%
Wrist 42.35% 18.49% 22.86% 10.23% 5.13% 0.78% 0.15%
Necklace 30.13% 17.57% 27.39% 15.00% 8.33% 1.33% 0.25%
Sleeve collar 97.44% 1.34% 0.85% 0.25% 0.11% 0.02% 0.00%
Collar pin 96.92% 1.60% 1.03% 0.30% 0.13% 0.02% 0.00%
Neck loop 99.69% 0.16% 0.10% 0.03% 0.01% 0.00% 0.00%
Hat pin 100.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00%
Glasses 99.99% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00%

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)

Model selection

Code
parameter <- "wearlikelihoodrating_social"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for wearlikelihoodrating_social
desc p.value
interaction of sex and wearing position 0.032
inclusion of sex 0.824
inclusion of sample location 0.129
inclusion of wearing position <0.001

The data supports Wearing Position, and its interaction with Sex. Sample is not supported by the data. The resulting model m4 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m4"
Model <- Models$model[Models$name == chosen_model][[1]]

summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: wearlikelihoodrating_social ~ wearing_position * demographics_sex +  
    (1 | record_id)
data:    data
subset:  demographics_sex != "Other"

 link  threshold nobs logLik   AIC     niter       max.grad cond.H 
 logit flexible  1152 -1810.43 3664.86 2853(19418) 1.96e-03 6.7e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 3.094    1.759   
Number of groups:  record_id 144 

Coefficients:
                                                   Estimate Std. Error z value
wearing_positionWrist                              -0.86456    0.31040  -2.785
wearing_positionNecklace                           -1.39967    0.30101  -4.650
wearing_positionSleeve collar                      -2.67891    0.30999  -8.642
wearing_positionCollar pin                         -2.25883    0.31001  -7.286
wearing_positionNeck loop                          -2.84337    0.31532  -9.017
wearing_positionHat pin                            -3.36561    0.31872 -10.560
wearing_positionGlasses                            -3.63839    0.33821 -10.758
demographics_sexMale                               -0.61176    0.41761  -1.465
wearing_positionWrist:demographics_sexMale          0.66775    0.43052   1.551
wearing_positionNecklace:demographics_sexMale       0.70436    0.42624   1.652
wearing_positionSleeve collar:demographics_sexMale  1.20096    0.43381   2.768
wearing_positionCollar pin:demographics_sexMale     0.02352    0.42906   0.055
wearing_positionNeck loop:demographics_sexMale      0.98738    0.43120   2.290
wearing_positionHat pin:demographics_sexMale        1.10079    0.43432   2.535
wearing_positionGlasses:demographics_sexMale        1.15958    0.45243   2.563
                                                   Pr(>|z|)    
wearing_positionWrist                               0.00535 ** 
wearing_positionNecklace                           3.32e-06 ***
wearing_positionSleeve collar                       < 2e-16 ***
wearing_positionCollar pin                         3.19e-13 ***
wearing_positionNeck loop                           < 2e-16 ***
wearing_positionHat pin                             < 2e-16 ***
wearing_positionGlasses                             < 2e-16 ***
demographics_sexMale                                0.14295    
wearing_positionWrist:demographics_sexMale          0.12090    
wearing_positionNecklace:demographics_sexMale       0.09844 .  
wearing_positionSleeve collar:demographics_sexMale  0.00563 ** 
wearing_positionCollar pin:demographics_sexMale     0.95629    
wearing_positionNeck loop:demographics_sexMale      0.02203 *  
wearing_positionHat pin:demographics_sexMale        0.01126 *  
wearing_positionGlasses:demographics_sexMale        0.01038 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -4.0338     0.3200 -12.607
Very Unlikely|Unlikely            -2.9256     0.3105  -9.423
Unlikely|Neutral                  -1.7688     0.3031  -5.836
Neutral|Likely                    -0.6461     0.2987  -2.163
Likely|Very likely                 1.2000     0.3017   3.978
Very likely|Extremely likely       3.2553     0.3490   9.329
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.08991832
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model, subset = "true")
sig_matrix <- 
  sig_matrix %>% mutate(term = str_remove(term, "demographics_sex"))


P5 <- 
  Significance_matrix(
sig_matrix %>% filter(!str_detect(term, "Male")), sex = TRUE
) +
geom_tile(
data = sig_matrix %>% filter(str_detect(term,"Male") & different) %>% 
        mutate(term = str_remove(term, ":Male")), 
fill = "red")

ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P5, width = 4, height = 4.7)

P5 

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
for(i in 1:144) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
abline(h=0, lty=2)

Model selection

Code
parameter <- "wearlikelihoodrating_exercise"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for wearlikelihoodrating_exercise
desc p.value
interaction of sex and wearing position 0.051
inclusion of sex 0.61
inclusion of sample location 0.506
inclusion of wearing position <0.001

The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: wearlikelihoodrating_exercise ~ wearing_position + (1 | record_id)
data:    data

 link  threshold nobs logLik   AIC     niter      max.grad cond.H 
 logit flexible  1160 -1915.38 3858.76 1662(8332) 1.14e-03 2.2e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 1.936    1.392   
Number of groups:  record_id 145 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
wearing_positionWrist          0.08459    0.21428   0.395    0.693    
wearing_positionNecklace      -1.97435    0.21822  -9.048  < 2e-16 ***
wearing_positionSleeve collar -0.12314    0.21650  -0.569    0.570    
wearing_positionCollar pin    -1.22176    0.21253  -5.749 8.99e-09 ***
wearing_positionNeck loop     -2.01988    0.21945  -9.204  < 2e-16 ***
wearing_positionHat pin       -2.52201    0.22273 -11.323  < 2e-16 ***
wearing_positionGlasses       -3.12963    0.23294 -13.435  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                                 Estimate Std. Error z value
Extremely unlikely|Very Unlikely  -3.6177     0.2210 -16.369
Very Unlikely|Unlikely            -2.7594     0.2108 -13.088
Unlikely|Neutral                  -1.6150     0.2003  -8.062
Neutral|Likely                    -0.7049     0.1947  -3.620
Likely|Very likely                 0.9452     0.1962   4.817
Very likely|Extremely likely       2.6629     0.2271  11.726
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.09162989
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model)

P6 <- Significance_matrix(sig_matrix)

ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P6, width = 4, height = 4.7)

P6

Predicted probabilities

Code
#if one wanted to print predictions for different subject percentiles
stDev <-  as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev, 
           "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_exercise
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 2.61% 3.34% 10.63% 16.48% 38.94% 21.46% 6.52%
Wrist 3.40% 4.26% 13.00% 18.62% 37.83% 17.83% 5.06%
Necklace 0.00% 0.01% 0.03% 0.05% 0.38% 2.12% 97.41%
Sleeve collar 1.78% 2.32% 7.74% 13.18% 38.45% 27.16% 9.36%
Collar pin 0.05% 0.07% 0.28% 0.59% 3.99% 17.65% 77.36%
Neck loop 0.00% 0.01% 0.02% 0.05% 0.33% 1.84% 97.75%
Hat pin 0.00% 0.00% 0.00% 0.01% 0.07% 0.38% 99.54%
Glasses 0.00% 0.00% 0.00% 0.00% 0.01% 0.06% 99.93%
Code
#average:
prob_table(Model$beta, 
           "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_exercise
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 2.61% 3.34% 10.63% 16.48% 38.94% 21.46% 6.52%
Wrist 2.41% 3.09% 9.95% 15.78% 39.05% 22.67% 7.06%
Necklace 16.20% 15.12% 27.56% 19.18% 16.82% 4.16% 0.96%
Sleeve collar 2.95% 3.74% 11.68% 17.49% 38.58% 19.76% 5.81%
Collar pin 8.35% 9.34% 22.60% 22.35% 27.08% 8.26% 2.01%
Neck loop 16.83% 15.48% 27.67% 18.85% 16.26% 3.99% 0.92%
Hat pin 25.05% 19.04% 27.14% 14.78% 10.95% 2.47% 0.56%
Glasses 38.04% 21.12% 22.82% 9.90% 6.46% 1.37% 0.30%
Code
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev, 
           "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_exercise
Wearing Position Extremely unlikely Very Unlikely Unlikely Neutral Likely Very likely Extremely likely
Chest Pin 2.61% 3.34% 10.63% 16.48% 38.94% 21.46% 6.52%
Wrist 2.01% 2.60% 8.57% 14.21% 38.88% 25.35% 8.37%
Necklace 93.53% 3.62% 1.92% 0.55% 0.30% 0.06% 0.01%
Sleeve collar 3.82% 4.75% 14.17% 19.50% 36.96% 16.29% 4.50%
Collar pin 56.80% 18.82% 15.07% 5.34% 3.18% 0.64% 0.14%
Neck loop 94.35% 3.17% 1.67% 0.48% 0.26% 0.05% 0.01%
Hat pin 98.80% 0.68% 0.35% 0.10% 0.05% 0.01% 0.00%
Glasses 99.83% 0.10% 0.05% 0.01% 0.01% 0.00% 0.00%

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)

Model selection

Code
parameter <- "wearduration"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for wearduration
desc p.value
interaction of sex and wearing position 0.005
inclusion of sex 0.852
inclusion of sample location 0.852
inclusion of wearing position <0.001

The data supports Wearing Position, and its interaction with Sex. Sample is not supported by the data. The resulting model m4 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m4"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: wearduration ~ wearing_position * demographics_sex + (1 | record_id)
data:    data
subset:  demographics_sex != "Other"

 link  threshold nobs logLik   AIC     niter       max.grad cond.H 
 logit flexible  1152 -1904.65 3853.30 3096(16129) 1.60e-03 7.1e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 2.521    1.588   
Number of groups:  record_id 144 

Coefficients:
                                                   Estimate Std. Error z value
wearing_positionWrist                               -0.7062     0.3044  -2.320
wearing_positionNecklace                            -1.4241     0.2999  -4.749
wearing_positionSleeve collar                       -2.3526     0.3103  -7.582
wearing_positionCollar pin                          -2.0286     0.3094  -6.557
wearing_positionNeck loop                           -3.2638     0.3185 -10.247
wearing_positionHat pin                             -4.0880     0.3266 -12.518
wearing_positionGlasses                             -3.8876     0.3300 -11.779
demographics_sexMale                                -0.8566     0.4047  -2.117
wearing_positionWrist:demographics_sexMale           0.5400     0.4336   1.245
wearing_positionNecklace:demographics_sexMale        0.3419     0.4268   0.801
wearing_positionSleeve collar:demographics_sexMale   1.1391     0.4337   2.627
wearing_positionCollar pin:demographics_sexMale      0.3680     0.4317   0.853
wearing_positionNeck loop:demographics_sexMale       1.3336     0.4350   3.066
wearing_positionHat pin:demographics_sexMale         1.5122     0.4376   3.456
wearing_positionGlasses:demographics_sexMale         1.1499     0.4450   2.584
                                                   Pr(>|z|)    
wearing_positionWrist                              0.020348 *  
wearing_positionNecklace                           2.04e-06 ***
wearing_positionSleeve collar                      3.39e-14 ***
wearing_positionCollar pin                         5.48e-11 ***
wearing_positionNeck loop                           < 2e-16 ***
wearing_positionHat pin                             < 2e-16 ***
wearing_positionGlasses                             < 2e-16 ***
demographics_sexMale                               0.034286 *  
wearing_positionWrist:demographics_sexMale         0.212968    
wearing_positionNecklace:demographics_sexMale      0.423084    
wearing_positionSleeve collar:demographics_sexMale 0.008626 ** 
wearing_positionCollar pin:demographics_sexMale    0.393928    
wearing_positionNeck loop:demographics_sexMale     0.002171 ** 
wearing_positionHat pin:demographics_sexMale       0.000549 ***
wearing_positionGlasses:demographics_sexMale       0.009758 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                              Estimate Std. Error z value
Not at all|Up to 1 hour        -4.6898     0.3164 -14.820
Up to 1 hour|1-3 hours         -3.6121     0.3061 -11.802
1-3 hours|3-6 hours            -2.2042     0.2950  -7.471
3-6 hours|6-8 hours            -1.0963     0.2886  -3.799
6-8 hours|8-12 hours           -0.1181     0.2855  -0.414
8-12 hours|Throughout the day   0.6863     0.2871   2.390
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.1005441
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model, subset = "true")
sig_matrix <- 
  sig_matrix %>% mutate(term = str_remove(term, "demographics_sex"))


P7 <- 
  Significance_matrix(
sig_matrix %>% filter(!str_detect(term, "Male")), sex = TRUE
) +
geom_tile(
data = sig_matrix %>% filter(str_detect(term,"Male") & different) %>% 
        mutate(term = str_remove(term, ":Male")), 
fill = "red")


ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P7, width = 4, height = 4.7)

P7

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
for(i in 1:144) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:144, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:144, labels=ord.re)
axis(2)
abline(h=0, lty=2)

Model selection

Code
parameter <- "restrict"
parameter_symbol <- sym(parameter)

#model generation
Models <- Infe(!!parameter_symbol, data_long)

#model comparisons
Model_Comparisons(Models) %>% 
  gt(caption = md(glue("Model comparisons for **{parameter}**"))) %>% 
  fmt_number(decimals = 3)
Model comparisons for restrict
desc p.value
interaction of sex and wearing position 0.3
inclusion of sex 0.918
inclusion of sample location 0.918
inclusion of wearing position <0.001

The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.

Significant differences

Code
#final model
chosen_model <- "m5"
Model <- Models$model[Models$name == chosen_model][[1]]
summary(Model)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
quadrature approximation with 10 quadrature points 

formula: restrict ~ wearing_position + (1 | record_id)
data:    data

 link  threshold nobs logLik   AIC     niter      max.grad cond.H 
 logit flexible  1160 -1764.02 3556.03 1587(4764) 9.38e-04 3.4e+02

Random effects:
 Groups    Name        Variance Std.Dev.
 record_id (Intercept) 1.201    1.096   
Number of groups:  record_id 145 

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
wearing_positionWrist          -1.3797     0.2457  -5.614 1.97e-08 ***
wearing_positionNecklace       -1.9653     0.2436  -8.069 7.11e-16 ***
wearing_positionSleeve collar  -1.4102     0.2523  -5.590 2.27e-08 ***
wearing_positionCollar pin     -1.9734     0.2489  -7.927 2.25e-15 ***
wearing_positionNeck loop      -2.3979     0.2478  -9.678  < 2e-16 ***
wearing_positionHat pin        -2.6988     0.2491 -10.836  < 2e-16 ***
wearing_positionGlasses        -3.1695     0.2561 -12.374  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                          Estimate Std. Error z value
Extremely|Significantly    -6.3272     0.3051 -20.735
Significantly|Quite a bit  -5.2148     0.2633 -19.807
Quite a bit|Moderately     -4.2970     0.2450 -17.540
Moderately|Somewhat        -3.4811     0.2335 -14.908
Somewhat|Slightly          -2.4176     0.2222 -10.882
Slightly|Not at all        -0.9715     0.2107  -4.612
Code
#McFadden
McFadden(Model, Models$model[[6]])
[1] 0.06010535
Code
#Significance matrix
sig_matrix <- data_significance_matrix(Model)

P8 <- Significance_matrix(sig_matrix)

ggsave(glue("output/04_unused/significance_matrix_{parameter}.pdf"), 
       P8, width = 4, height = 4.7)

P8

Predicted probabilities

Code
#if one wanted to print predictions for different subject percentiles
stDev <-  as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2
# 5th percentile:
prob_table(Model$beta * qnorm(0.05) * stDev, 
           "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for restrict
Wearing Position Extremely Significantly Quite a bit Moderately Somewhat Slightly Not at all
Chest Pin 0.18% 0.36% 0.80% 1.64% 5.20% 19.27% 72.54%
Wrist 0.01% 0.02% 0.05% 0.11% 0.38% 1.84% 97.58%
Necklace 0.00% 0.01% 0.02% 0.04% 0.12% 0.59% 99.23%
Sleeve collar 0.01% 0.02% 0.05% 0.11% 0.36% 1.74% 97.72%
Collar pin 0.00% 0.01% 0.02% 0.03% 0.12% 0.58% 99.24%
Neck loop 0.00% 0.00% 0.01% 0.02% 0.05% 0.25% 99.67%
Hat pin 0.00% 0.00% 0.00% 0.01% 0.03% 0.14% 99.82%
Glasses 0.00% 0.00% 0.00% 0.00% 0.01% 0.06% 99.93%
Code
#average:
prob_table(Model$beta, 
           "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for restrict
Wearing Position Extremely Significantly Quite a bit Moderately Somewhat Slightly Not at all
Chest Pin 0.18% 0.36% 0.80% 1.64% 5.20% 19.27% 72.54%
Wrist 0.71% 1.41% 3.02% 5.77% 15.26% 33.91% 39.94%
Necklace 1.26% 2.48% 5.12% 9.16% 20.87% 34.10% 27.02%
Sleeve collar 0.73% 1.45% 3.10% 5.91% 15.55% 34.05% 39.21%
Collar pin 1.27% 2.49% 5.15% 9.21% 20.95% 34.07% 26.86%
Neck loop 1.93% 3.71% 7.38% 12.27% 24.22% 31.13% 19.37%
Hat pin 2.59% 4.89% 9.35% 14.56% 25.60% 27.92% 15.09%
Glasses 4.08% 7.37% 13.01% 17.81% 25.69% 22.05% 9.99%
Code
# 95th percentile:
prob_table(Model$beta * qnorm(0.95) * stDev, 
           "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for restrict
Wearing Position Extremely Significantly Quite a bit Moderately Somewhat Slightly Not at all
Chest Pin 0.18% 0.36% 0.80% 1.64% 5.20% 19.27% 72.54%
Wrist 2.65% 5.00% 9.54% 14.76% 25.68% 27.61% 14.76%
Necklace 7.98% 12.89% 18.90% 20.12% 21.33% 13.62% 5.16%
Sleeve collar 2.81% 5.28% 9.97% 15.21% 25.82% 26.89% 14.02%
Collar pin 8.10% 13.04% 19.02% 20.12% 21.19% 13.45% 5.09%
Neck loop 16.93% 21.34% 22.55% 17.01% 13.22% 6.69% 2.26%
Hat pin 26.96% 25.93% 20.87% 12.64% 8.44% 3.89% 1.26%
Glasses 48.33% 25.66% 13.69% 6.46% 3.75% 1.60% 0.50%

Model diagnostics

Code
#random effect distribution with random effects
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
for(i in 1:145) segments(i, ci[i,1], i, ci[i,2])
abline(h=0, lty=2)

Code
#without confidence intervals
ci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o% c(-1,1)
ord.re <- order(Model$ranef)
ci <- ci[order(Model$ranef),]
plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")
axis(1, at=1:145, labels=ord.re)
axis(2)
abline(h=0, lty=2)

Summary Graph

Code
a <- theme(plot.margin = margin(5, 10,5,5,"pt"))
b <- theme(plot.margin = margin(5, 5,5,10,"pt"))

(P1 + a + P2 + b) / (P3 + a + P4 + b) / (P5 + a + P6 + b) / (P7 + a + P8 + b) + 
  plot_annotation(tag_levels = "A") & theme(plot.tag = element_text(size = 20))

Code
ggsave("output/01_figures/figure2_summary_graph_differences.pdf", width = 10, height = 20)
ggsave("output/01_figures/figure2_summary_graph_differences.jpeg", width = 10, height = 20, dpi = 300)